#include <stdlib.h>
#include <fftw3.h>

//傅里叶变换
void FFT2(double *temp , int Im_W , int Im_H)
{
	int i,j;
	fftw_complex *in, *out;
    fftw_plan p;
    //分配存储空间
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * Im_H*Im_W);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * Im_H*Im_W);
    //设置变换计划
    p = fftw_plan_dft_2d(Im_H, Im_W, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    for (i=0; i<Im_H; i++){
        for (j=0; j<Im_W; j++){
            in[i*Im_W+j][0] = *(temp+i*Im_W*2+j*2+0);
            in[i*Im_W+j][1] = *(temp+i*Im_W*2+j*2+1);
        }
    }
    fftw_execute(p);// 执行变换

	for (i=0; i<Im_H; i++){
        for (j=0; j<Im_W; j++){
				*(temp+i*Im_W*2+j*2+0) = out[i*Im_W+j][0];
				*(temp+i*Im_W*2+j*2+1) = out[i*Im_W+j][1];
			}
    }
    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);
}
//逆傅里叶变换
void iFFT2(double *temp , int Im_W , int Im_H)
{
	int i,j;
	fftw_complex *in, *out;
    fftw_plan p;
    //分配存储空间
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * Im_H*Im_W);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * Im_H*Im_W);
    //设置变换计划
    p = fftw_plan_dft_2d(Im_H, Im_W, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
    for (i=0; i<Im_H; i++){
        for (j=0; j<Im_W; j++){
            in[i*Im_W+j][0] = *(temp+i*Im_W*2+j*2+0);
            in[i*Im_W+j][1] = *(temp+i*Im_W*2+j*2+1);
        }
    }
    fftw_execute(p);//执行变换

	for (i=0; i<Im_H; i++){
        for (j=0; j<Im_W; j++){
				*(temp+i*Im_W*2+j*2+0) = out[i*Im_W+j][0]/Im_H/Im_W;
				*(temp+i*Im_W*2+j*2+1) = out[i*Im_W+j][1]/Im_H/Im_W;
			}
    }
    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);
}
//移动谱线级次
void fftshift(double *temp , int Im_W , int Im_H)
{
	int i,j,dl, a = Im_H/2 , b = Im_W/2;
	double *temp_s = (double*)malloc(sizeof(double)*Im_W*Im_H*2);
	for(i=0 ; i<Im_H ; i++){
		dl = (i+a)%Im_H;
		for(j=0 ; j<Im_W ; j++){
			*(temp_s + dl*Im_W*2 + j*2 + 0) = *(temp + i*Im_W*2 + j*2 +0);
			*(temp_s + dl*Im_W*2 + j*2 + 1) = *(temp + i*Im_W*2 + j*2 +1);
		}
	}
	for(j=0 ; j<Im_W ; j++){
		dl = (j+b)%Im_W;
		for(i=0 ; i< Im_H ; i++){
			*(temp + i*Im_W*2 + dl*2 +0) = *(temp_s + i*Im_W*2 + j*2 +0);
			*(temp + i*Im_W*2 + dl*2 +1) = *(temp_s + i*Im_W*2 + j*2 +1);
		}
	}

	free(temp_s);
}
